This tutorial illustrates how to use a database of strong ground motions to expected ground motion values from GMPEs with observations.
A database (of metadata) is stored as a Python "Pickle" file. To load the database, do the following:
In [ ]:
%load_ext autoreload
%autoreload 2
import warnings; warnings.filterwarnings("ignore")
In [ ]:
%matplotlib inline
import numpy as np
import cPickle
input_database = "data/demo_records/metadata.pkl"
db1 = cPickle.load(open(input_database, "r"))
print "Number of records: %g" % db1.number_records()
In the following example we select a subset of records from a database and compare these observed ground motions against four relevent GMPEs:
1.) Boore & Atkinson (2008) 2.) Akkar & Bommer (2010) 3.) Akkar & Cagnan (2010) 4.) Akkar et al (2014) - Joyner-Boore coefficients
Four intensity measures will also be considered: PGA, Sa (0.2), Sa(1.0), Sa(2.0)
In [ ]:
gmpe_list = ["BooreAtkinson2008",
"AkkarBommer2010",
"AkkarCagnan2010",
"AkkarEtAlRjb2014"]
imt_list = ["PGA", "SA(0.2)", "SA(1.0)", "SA(2.0)"]
In [ ]:
# Get residuals
from smtk.residuals.gmpe_residuals import Residuals, Likelihood, LLH, EDR
import smtk.residuals.residual_plotter as rspl
# Create an instance of the residual calculator
resid1 = Residuals(gmpe_list, imt_list)
resid1.get_residuals(db1)
In [ ]:
residuals_file = "data/demo_residuals1.csv"
resid1.pretty_print(residuals_file, sep=",")
In [ ]:
rspl.ResidualPlot(resid1, "BooreAtkinson2008", "PGA")
In [ ]:
rspl.ResidualPlot(resid1, "AkkarEtAlRjb2014", "PGA")
In [ ]:
lh_analysis = Likelihood(gmpe_list, imt_list)
lh_analysis.get_residuals(db1)
In [ ]:
rspl.LikelihoodPlot(lh_analysis, "BooreAtkinson2008", "PGA")
In [ ]:
rspl.LikelihoodPlot(lh_analysis, "AkkarEtAlRjb2014", "PGA")
In [ ]:
rspl.LikelihoodPlot(lh_analysis, "AkkarEtAlRjb2014", "SA(1.0)")
In [ ]:
rspl.ResidualWithDistance(resid1, "AkkarEtAlRjb2014", "PGA", plot_type="linear", figure_size=(7, 7))
In [ ]:
rspl.ResidualWithDistance(resid1, "BooreAtkinson2008", "PGA", plot_type="linear", figure_size=(7, 7))
In [ ]:
rspl.ResidualWithMagnitude(resid1, "AkkarEtAlRjb2014", "PGA", figure_size=(7,7))
In [ ]:
rspl.ResidualWithMagnitude(resid1, "BooreAtkinson2008", "PGA", figure_size=(7,7))
In [ ]:
rspl.ResidualWithVs30(resid1, "AkkarEtAlRjb2014", "PGA", figure_size=(7,7))
In [ ]:
rspl.ResidualWithVs30(resid1, "BooreAtkinson2008", "PGA", figure_size=(7,7))
In [ ]:
rspl.ResidualWithDepth(resid1, "AkkarEtAlRjb2014", "PGA", figure_size=(7,7))
In [ ]:
rspl.ResidualWithDepth(resid1, "BooreAtkinson2008", "PGA", figure_size=(7,7))
In [ ]:
# Create an instance of the residual calculator
llh_1 = LLH(gmpe_list, imt_list)
llh_1.get_residuals(db1)
In [ ]:
# Run LLH analysis
llh_values, weights = llh_1.get_loglikelihood_values(imt_list)
for gmpe in gmpe_list:
print "LLH(%s) = %12.6f Weight = %12.6f" % (gmpe, llh_values[gmpe]["All"], weights[gmpe])
In [ ]:
# EDR analysis
edr1 = EDR(gmpe_list, imt_list)
edr1.get_residuals(db1)
In [ ]:
edr_values = edr1.get_edr_values()
for gmpe in gmpe_list:
print "EDR(%s) = %8.4f (sqrt(kappa) = %8.4f)" %(gmpe,
edr_values[gmpe]["EDR"],
edr_values[gmpe]["sqrt Kappa"])
This illustrates the application of a set of site specific analyses to determine the single-station phi. Implementation is based on the description and formulae found in Rodgriguez-Marek et al. (2011)
In [ ]:
from smtk.strong_motion_selector import SMRecordSelector, rank_sites_by_record_count
import smtk.residuals.gmpe_residuals as res
import smtk.residuals.residual_plotter as rspl
In [ ]:
# Returns a list of the stations with the most records (in descending order)
# cutting off sites with less than 30 records
top_sites = rank_sites_by_record_count(db1, threshold=5)
for site_id in top_sites.keys():
print "Site ID: %s Name: %s Number of Records: %s" %(site_id, top_sites[site_id]["Name"], top_sites[site_id]["Count"])
In [ ]:
# Consider the AkkarEtAlRjb2014 GMPE with PGA, Sa (0.2), Sa(1.0)
gmpe_list = ["AkkarEtAlRjb2014", "BooreEtAl2014"]
imt_list = ["PGA", "SA(0.2)", "SA(1.0)"]
# Get the residuals for the site
ssa1 = res.SingleStationAnalysis(top_sites.keys(), gmpe_list, imt_list)
ssa1.get_site_residuals(db1)
In [ ]:
# Plots the GMPE Residuals for Each Site
rspl.ResidualWithSite(ssa1, "AkkarEtAlRjb2014", "PGA", figure_size=(7,7))
In [ ]:
rspl.ResidualWithSite(ssa1, "AkkarEtAlRjb2014", "SA(1.0)", figure_size=(7,7))
In [ ]:
# Plots the Intra-event analysis for each site
rspl.IntraEventResidualWithSite(ssa1, "AkkarEtAlRjb2014", "PGA", figure_size=(7,7))
In [ ]:
rspl.IntraEventResidualWithSite(ssa1, "BooreEtAl2014", "PGA", figure_size=(7,7))
In [ ]:
# Plots the Intra-event analysis for each site
rspl.IntraEventResidualWithSite(ssa1, "AkkarEtAlRjb2014", "SA(1.0)", figure_size=(7,7))
In [ ]:
rspl.IntraEventResidualWithSite(ssa1, "BooreEtAl2014", "SA(1.0)", figure_size=(7,7))
In [ ]:
for iloc, output in enumerate(ssa1.site_residuals):
print "Site ID = %s" % (top_sites.keys()[iloc])
for gmpe in output.site_analysis.keys():
print "GMPE = %s" % gmpe
for imt_type in output.site_analysis[gmpe]:
site_results = output.site_analysis[gmpe][imt_type]
print "%10s: dS2ss = %10.8f phi_ss,s = %10.8f" % (imt_type, site_results["dS2ss"], site_results["phi_ss,s"])
print "======================================================="